#########################################################################################
#############################    Preparing Simulation     ###############################
#########################################################################################
#########################################################################################

library(pacman)
p_load(tidyverse, msm, sandwich, readxl)
set.seed(12345)
rm(list = ls())

#setwd('`~/Dropbox/StrengthInNumbersReplicationPackage/replicable')

#########################################################################################


load(file = "simulations/simulations_output_1.Rda")


modelInputsParams <- matrix(1:12, ncol = 3)


# set functions
getX <- function(y, model) {
  b <- coef(model)[1]
  c <- coef(model)[2]
  d <- coef(model)[3]
  x <- -(sqrt(-4*b*d+ c^2 + 4*d*y) + c)/(2*d)
  return(x)
}

self_voting <- read.csv("figure_5/figure_5_input.csv")

femMajValues  <- c(self_voting$X_val[18]/self_voting$X_exp[18])
maleMajValues <- c(self_voting$X_val[13]/self_voting$X_exp[13])

n_group_femmaj <- self_voting$group_N[8]
n_group_malmaj <- self_voting$group_N[3]


femMajDiscrim <- getX(femMajValues, femaleMajorityTrend)
maleMajDiscrim <- getX(maleMajValues, maleMajorityTrend)


pointsFrame <- as.data.frame(cbind(femMajDiscrim, femMajValues, maleMajDiscrim, maleMajValues))
outcomes <- c("Most influential self votes per woman\n(experimental mean)\n")
pointsFrame <- cbind(pointsFrame, outcomes)



estmean <- coef(femaleMajorityTrend_inverted)

estvar <- vcovCL(femaleMajorityTrend_inverted, cluster = ~`Performance penalty for women on themselves (SD)`) 
discrim1<- pointsFrame$femMajDiscrim[1]
pointsFrame$femMajDiscrimSE <- c(deltamethod (~ x1 + x2 * discrim1 + x3*discrim1*discrim1,
                                              estmean, estvar))
# now for men
estmean <- coef(maleMajorityTrend)
estvar <- vcovCL(maleMajorityTrend_inverted, cluster = ~`Performance penalty for women on themselves (SD)`)
discrim1<- pointsFrame$maleMajDiscrim[1]
pointsFrame$maleMajDiscrimSE <- c(deltamethod (~ x1 + x2 * discrim1 + x3*discrim1*discrim1,
                                               estmean, estvar))
pointsFrame

modelInputsParams[1,2]<-n_group_malmaj
modelInputsParams[1,3]<-n_group_femmaj
modelInputsParams[2,2]<-maleMajValues
modelInputsParams[2,3]<-femMajValues
modelInputsParams[3,2]<-maleMajDiscrim
modelInputsParams[3,3]<-femMajDiscrim
modelInputsParams[4,2]<-pointsFrame$maleMajDiscrimSE
modelInputsParams[4,3]<-pointsFrame$femMajDiscrimSE
modelInputsParams <- as.data.frame(modelInputsParams) 


modelInputsParams$V1[1] <- "Number of female-majority groups"
modelInputsParams$V1[2] <- "Women's rate of self-voting (normalized)"
modelInputsParams$V1[3] <- "Women's implied self-penalty (in SD)"
modelInputsParams$V1[4] <- "(se)"

colnames(modelInputsParams) <- c("Measure", "Male-majority groups", "Female-majority groups")

save(modelInputsParams, pointsFrame, file = "simulations/simulations_input_params.Rda")

